There are two main types of spatial data: vector and raster data.
Vector: Vectors (also called shapefiles) consist of points, lines and polygons. These shapes are attached to a dataframe, where each row corresponds to a different spatial element.
Raster: Rasters are spatially-referenced grids where each cell has one value.
Shapefile
Raster
Here, we'll use the following packages
# Packages and Filepaths ------------------------------------------------------- library(sp) library(rgdal) library(rgeos) library(raster) hh_data <- read.csv(file.path(finalData,"HH_data.csv")) str(hh_data)
## 'data.frame': 422 obs. of 5 variables: ## $ X : int 1 3 4 5 6 7 8 10 11 12 ... ## $ expend_food_yearly: num 250452 420029 484468 326109 131487 ... ## $ food_security : Factor w/ 3 levels "Little Hunger",..: 1 2 1 1 3 3 1 2 3 1 ... ## $ longitude_scramble: num 30.4 30.4 30.4 30.4 30.4 ... ## $ latitude_scramble : num -1.97 -2 -1.99 -1.98 -1.95 ...
library(ggplot2)
ggplot() +
geom_point(data=hh_data,
aes(x=longitude_scramble, y=latitude_scramble, color=food_security),
size=.7)
hh_map1 <- ggplot() +
# Points
geom_point(data=hh_data,
aes(x=longitude_scramble, y=latitude_scramble, color=food_security),
size=.7) +
# Other elements to improve map
coord_quickmap() + # make sure the map doesn't look distorted. Can also use
# coord_map(), but sometimes makes the process slow
theme_void() +
scale_color_manual(values=c("green", "orange", "red")) +
labs(title="Household Food Security", color="Food Security") +
theme(plot.title = element_text(hjust = 0.5, face="bold")) # center and bold title
hh_map1
library(ggmap)
basemap <- get_map(location = c(lon=mean(hh_data$longitude_scramble),
lat=mean(hh_data$latitude_scramble)),
zoom=10,
maptype="roadmap") # roadmap, satellite, etc. See help(get_map)
hh_map2 <- ggmap(basemap) +
geom_point(data=hh_data,
aes(x=longitude_scramble,
y=latitude_scramble,
color=food_security),
size=.7) +
coord_quickmap() + # make sure the map doesn't look distorted. Can also use
# coord_map(), but sometimes makes the process slow
theme_void() +
labs(title="Household Food Security", color="Food Security") +
theme(plot.title = element_text(hjust = 0.5, face="bold")) +
scale_color_manual(values=c("green", "orange", "red"))
hh_map2
hh_map2 +
scale_x_continuous(limits = c(min(hh_data$longitude_scramble),
max(hh_data$longitude_scramble)),
expand = c(.03, .03)) +
scale_y_continuous(limits = c(min(hh_data$latitude_scramble),
max(hh_data$latitude_scramble)),
expand = c(.03, .03))
So far, we've just been working with a dataframe and using the longitude and latitude variables. However, for many other uses of spatial data we need to convert the dataframe into a spatial dataframe.
coordinates transforms a dataframe into a spatial dataframe. The syntax is: coordinates(DATA FRAME) <- ~LONGITUDE+LATITUDE, where LONGITUDE and LATITUDE are the names of variables in DATA FRAME.crs.hh_data_sdf <- hh_data
coordinates(hh_data_sdf) <- ~longitude_scramble+latitude_scramble
crs(hh_data_sdf) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
Here, the variables are separated from the coordinates. Specifically, a spatial dataframe is a list, where one element of the list is a dataframe of variables and other is a dataframe of coordinates.
hh_data_sdf
## class : SpatialPointsDataFrame ## features : 422 ## extent : 30.26138, 30.67035, -2.094064, -1.934281 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ## variables : 3 ## names : X, expend_food_yearly, food_security ## min values : 1, 0.000, Little Hunger ## max values : 553, 3034121.500, Severe Hunger
Access a dataframe containing the variables (except the coordinates)
hh_data_sdf@data
Access the coordinates of the points
hh_data_sdf@coords
We can still access variables like a normal dataframe
hh_data_sdf@data$food_security hh_data_sdf$food_security # short-cut
We can quickly plot spatial dataframes.
plot(hh_data_sdf)
See rgeos for a list of common functions to apply on vector data.
library(rgeos) hh_data_sdf_buff <- gBuffer(hh_data_sdf, width= 0.5/111, byid=T) #buffer by about 10km plot(hh_data_sdf_buff)
head(hh_data_sdf@coords)
## longitude_scramble latitude_scramble ## 1 30.41310 -1.965660 ## 2 30.42626 -2.003319 ## 3 30.42944 -1.986508 ## 4 30.42331 -1.976507 ## 5 30.42285 -1.952778 ## 6 30.41069 -2.006962
hh_data_sdf_newproj <- spTransform(hh_data_sdf, CRS("+init=epsg:3857"))
head(hh_data_sdf_newproj@coords)
## longitude_scramble latitude_scramble ## 1 3385571 -218859.2 ## 2 3387035 -223053.9 ## 3 3387390 -221181.4 ## 4 3386707 -220067.4 ## 5 3386657 -217424.3 ## 6 3385302 -223459.7
Note: See appendix for ways to calculate distances in a more accurate way (e.g., by using an equal area projection or by taking into account the curvator of the earth)
dist_matrix <- gDistance(hh_data_sdf_newproj, byid=T) dist_matrix[1:7,1:7]
## 1 2 3 4 5 6 7 ## 1 0.000 4442.9586 2950.070 1658.693 1799.332 4608.262 4019.3068 ## 2 4442.959 0.0000 1905.795 3004.456 5642.253 1779.875 850.3531 ## 3 2950.070 1905.7947 0.000 1306.708 3828.044 3090.272 1193.2667 ## 4 1658.693 3004.4565 1306.708 0.000 2643.513 3671.725 2442.9956 ## 5 1799.332 5642.2528 3828.044 2643.513 0.000 6185.365 5019.0051 ## 6 4608.262 1779.8748 3090.272 3671.725 6185.365 0.000 2486.2690 ## 7 4019.307 850.3531 1193.267 2442.996 5019.005 2486.269 0.0000
So far, we've been making static maps. But those are boring. Let's make an interactive map using Leaflet. There's a bunch of different basemaps to choose from; from the list here, enter the name of the basemap in quotes in the addProviderTiles function.
library(leaflet)
library(dplyr)
imap_1 <- leaflet() %>%
addProviderTiles("OpenStreetMap") %>%
addCircles(data=hh_data_sdf)
NOTE: leaflet assumes you are using the following coordinate reference system:
"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
imap_1
We can save our interactive map outside of R. Here, we use saveWidget from the htmlwidgets package. Doubleclicking the html file will open the file in any browser.
library(htmlwidgets) saveWidget(imap_1, file=file.path(Output,"interactive_map.html"))
# Define Palette
pal <- colorFactor(
palette = c("Green","Yellow","Red"),
domain = hh_data_sdf$food_security)
imap_2 <- leaflet() %>%
addProviderTiles("Stamen.Terrain") %>%
addCircleMarkers(data=hh_data_sdf,
radius=3,
fillOpacity=1,
color=~pal(food_security),
stroke=F, # remove outer-circle around circle
popup=~food_security) %>% # variable to display when click feature
# Add legend for points layer
addLegend(pal = pal,
values = hh_data_sdf$food_security,
title = "Food Security")
imap_2
# Load and Prep Plot-Level Data ------------------------------------------------ setwd(file.path(finalData,"Shapefiles")) # set filepath to folder with shapefile ag_fields <- readOGR(dsn=".", layer="allsitessubset", verbose = FALSE) # Plot File Projection crs(ag_fields)
## CRS arguments: ## +proj=tmerc +lat_0=0 +lon_0=30 +k=0.9999 +x_0=500000 +y_0=5000000 ## +ellps=GRS80 +units=m +no_defs
# HH Survey Location Projection crs(hh_data_sdf)
## CRS arguments: ## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
# Reproject plot data to have same projection as HH survey data hh_dd_proj <- as.character(crs(hh_data_sdf)) ag_fields <- spTransform(ag_fields, CRS(hh_dd_proj))
plot(ag_fields)
ag_fields
## class : SpatialPolygonsDataFrame ## features : 19 ## extent : 30.25664, 30.67621, -2.095173, -1.930931 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ## variables : 1 ## names : site ## min values : Kayonza15 ## max values : Rwamagana35
ag_fields@polygons Each polyon is defined by a matrix of vertices
head(ag_fields)
## site ## 0 Rwamagana33 ## 1 Rwamagana33 ## 2 Kayonza15 ## 3 Kayonza15 ## 4 Rwamagana34 ## 5 Rwamagana34
You'll notice that multiple features have the same name. Here, let's simplify the polygon so that each site name represents one feature.

To aggregate (or dissolve) the polygon by the variable site, use the aggregate function from the raster package.
ag_fields <- aggregate(ag_fields, by = "site")
plot(ag_fields)
In the previous slide, you'll notice there's a weird tic mark in one of the polygons. To get rid of this, let's buffer the polygon by a very small amount – which will help cover the hole in the field.
ag_fields <- gBuffer(ag_fields, width=0.000001/111, byid=T) plot(ag_fields)
Before we make a map, let's merge in plot some data for each agriculture field. I've already prepared this data, but it includes average food securtiy and average food expenditure per agriculture field based on the household data.
ag_fields_data <- read.csv(file.path(finalData, "plot_data.csv")) ag_fields <- merge(ag_fields, ag_fields_data, by="site") summary(ag_fields@data)
## site expend_food_yearly_mean food_security_mean ## Kayonza15 :1 Min. :189101 Min. :1.468 ## Kayonza4 :1 1st Qu.:262529 1st Qu.:1.599 ## Rwamagana2 :1 Median :321472 Median :1.736 ## Rwamagana33:1 Mean :304092 Mean :1.726 ## Rwamagana34:1 3rd Qu.:354403 3rd Qu.:1.898 ## Rwamagana35:1 Max. :384243 Max. :1.915
ggplot can't interpret spatial dataframes. Consequently, we need to transform the data into a format that ggplot can understand. Here, we make a dataframe where each vertice of the polygon is an observation.
library(broom) ag_fields_df <- tidy(ag_fields, region="site") head(ag_fields_df)
## long lat order hole piece group id ## 1 30.58359 -1.936452 1 FALSE 1 Kayonza15.1 Kayonza15 ## 2 30.58364 -1.936513 2 FALSE 1 Kayonza15.1 Kayonza15 ## 3 30.58364 -1.936496 3 FALSE 1 Kayonza15.1 Kayonza15 ## 4 30.58366 -1.936071 4 FALSE 1 Kayonza15.1 Kayonza15 ## 5 30.58366 -1.936071 5 FALSE 1 Kayonza15.1 Kayonza15 ## 6 30.58366 -1.936032 6 FALSE 1 Kayonza15.1 Kayonza15
The resulting dataframe from tidy:
id, which is taken from the the variable in the region argument.ag_fields) into this dataframe using the site variable from ag_fields and the id variable from ag_fields_df.ag_fields_df <- merge(ag_fields_df, ag_fields, by.x="id", by.y="site") head(ag_fields_df)
## id long lat order hole piece group ## 1 Kayonza15 30.58359 -1.936452 1 FALSE 1 Kayonza15.1 ## 2 Kayonza15 30.58364 -1.936513 2 FALSE 1 Kayonza15.1 ## 3 Kayonza15 30.58364 -1.936496 3 FALSE 1 Kayonza15.1 ## 4 Kayonza15 30.58366 -1.936071 4 FALSE 1 Kayonza15.1 ## 5 Kayonza15 30.58366 -1.936071 5 FALSE 1 Kayonza15.1 ## 6 Kayonza15 30.58366 -1.936032 6 FALSE 1 Kayonza15.1 ## expend_food_yearly_mean food_security_mean ## 1 189101.4 1.914634 ## 2 189101.4 1.914634 ## 3 189101.4 1.914634 ## 4 189101.4 1.914634 ## 5 189101.4 1.914634 ## 6 189101.4 1.914634
ggplot() +
geom_polygon(data = ag_fields_df, aes(x = long, y=lat, group=group,
fill=expend_food_yearly_mean))
ag_map1 <- ggplot() +
geom_polygon(data = ag_fields_df, aes(x = long, y=lat, group=group,
fill=expend_food_yearly_mean),
color="black", size=0.25) + # Borders of polygons
coord_quickmap() + # Make sure map doesn't look distorted
theme_void() + # Set theme of plot
scale_fill_gradient(low = "orangered1", high = "green1") + # Color Gradient
labs(fill="Food\nExpenditure", title="Annual Food Expenditure") + # Labels
theme(plot.title = element_text(hjust = 0.5, face="bold")) # Center title
ag_map1
library(RColorBrewer)
ag_map2 <- ggplot() +
geom_polygon(data = ag_fields_df, aes(x = long, y=lat, group=group,
fill=expend_food_yearly_mean),
color="black", size=0.25) + # Borders of polygons
coord_quickmap() +
theme_void() +
scale_fill_distiller(palette = "Spectral", direction = 1) + # Color Gradient
labs(fill="Food\nExpenditure", title="Annual Food Expenditure") +
theme(plot.title = element_text(hjust = 0.5, face="bold"))
ag_map2
scale_*_gradient(low, high): [For continuous variable] Manually define low/high colorsscale_*_gradientn(colors): [For continuous variable] Use a defined list of colors (e.g., c("purple","blue","yellow","white"))scale_*_manual(colors): [For discrete variable] Use a defined list of colors, where the list is equal to the number of unique observations in the discrete variable.scale_*_distiller(palette): [For contiuous variables]scale_*_brewer(palette): [For discrete variables]Where * is either color or fill. Discrete variables = factor variables.
Now, let's add text to the map showing regions. To do this, we need to create a dataframe that includes the lat/lon of the center of each plot and a variable for the site name. To do this, we use:
gCentroid() which returns a shapefile of the center of each regioncoordinates which returns a matrix of the coordinates of a shapefile, where the first column is longitude and the second column is latitude.# Create dataframe of center of each plot with site name as variable
ag_fields_center <- gCentroid(ag_fields, byid=T)
ag_fields_center <- as.data.frame(coordinates(ag_fields_center))
ag_fields_center$site <- ag_fields$site
names(ag_fields_center) <- c("longitude","latitude","site")
Now, lets use geom_text to add text to our map.
ag_map2 + geom_text(data=ag_fields_center, aes(label=site, x=longitude, y=latitude),
check_overlap = TRUE) # makes sure text doesn't overlap
If you have lots of cluttered text, use geom_text_repel from ggrepel.
library(ggrepel) ag_map2 + geom_text_repel(data=ag_fields_center, aes(label=site, x=longitude, y=latitude))
ag_map2 <- ggplot() +
geom_polygon(data = ag_fields_df, aes(x = long, y=lat, group=group,
fill=expend_food_yearly_mean),
color="black", size=0.25) + # Borders of polygons
geom_point(data=hh_data,
aes(x=longitude_scramble,
y=latitude_scramble,
color="HH Location"), # Name an aesthetic want you want to
# appear on legend
size=.1, alpha=.6) +
scale_color_manual(values="black") + # Manually define color of points
coord_quickmap() +
theme_void() +
scale_fill_gradient(low = "orangered1", high = "green1") +
labs(fill="Food\nExpenditure", title="Annual Food Expenditure",
color="") + # Setting color="" makes the title above the legend item blank
theme(plot.title = element_text(hjust = 0.5, face="bold"))
ag_map2
leaflet() %>%
addProviderTiles("OpenStreetMap") %>%
addPolygons(data=ag_fields)
pal_foodexpend <- colorNumeric("RdYlGn", ag_fields$expend_food_yearly_mean)
imap_3 <- leaflet() %>%
addProviderTiles("OpenStreetMap") %>%
addPolygons(data=ag_fields,
fillColor = ~pal_foodexpend(expend_food_yearly_mean),
fillOpacity = 0.6,
color="black", # color of line around polygons
weight=1, # width of line around polygons
popup=~site)
imap_3
Now, let's make a map using administrative-level data. The get_data function from the raster package allows us to download data from GADM (the Database of Global Administrative Areas). Let's download the third administrative division for Rwanda.
rwa_adm <- getData('GADM', country='RWA', level=3)
rwa_adm
## class : SpatialPolygonsDataFrame ## features : 422 ## extent : 28.86171, 30.89907, -2.839973, -1.04745 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ## variables : 16 ## names : OBJECTID, ID_0, ISO, NAME_0, ID_1, NAME_1, ID_2, NAME_2, ID_3, NAME_3, CCN_3, CCA_3, TYPE_3, ENGTYPE_3, NL_NAME_3, ... ## min values : 1, 189, RWA, Rwanda, 1, Amajyaruguru, 1, Bugesera, 1, Base, 0, 0, Sector, Sector, , ... ## max values : 422, 189, RWA, Rwanda, 5, Umujyi wa Kigali, 30, Rwamagana, 422, Zaza, 5715, 5715, Water body, Water body, , ...
head(rwa_adm)
## OBJECTID ID_0 ISO NAME_0 ID_1 NAME_1 ID_2 NAME_2 ID_3 NAME_3 ## 1 1 189 RWA Rwanda 1 Amajyaruguru 1 Burera 1 Bungwe ## 2 2 189 RWA Rwanda 1 Amajyaruguru 1 Burera 2 Butaro ## 3 3 189 RWA Rwanda 1 Amajyaruguru 1 Burera 3 Cyanika ## 4 4 189 RWA Rwanda 1 Amajyaruguru 1 Burera 4 Cyeru ## 5 5 189 RWA Rwanda 1 Amajyaruguru 1 Burera 5 Gahunga ## 6 6 189 RWA Rwanda 1 Amajyaruguru 1 Burera 6 Gatebe ## CCN_3 CCA_3 TYPE_3 ENGTYPE_3 NL_NAME_3 VARNAME_3 ## 1 4401 4401 Sector Sector ## 2 4402 4402 Sector Sector ## 3 4403 4403 Sector Sector ## 4 4404 4404 Sector Sector ## 5 4405 4405 Sector Sector ## 6 4406 4406 Sector Sector
I've preppared a dataset of average vegetation levels per administrative unit in Rwanda. Here, 0 means no vegetation and 1 means significant vegetation.
rwa_veg <- read.csv(file.path(finalData, "rwanda_vegetation.csv")) head(rwa_veg)
## X OBJECTID vegetation ## 1 1 1 0.542 ## 2 2 2 0.602 ## 3 3 3 0.677 ## 4 4 4 0.604 ## 5 5 5 0.591 ## 6 6 6 0.546
Make an administrative-level map of vegetation levels. Here's a broad outline of steps you'll need to implement:
tol make the polygon more simplified. I played around with values and found 0.003 to be adequate. Note that large values of tol might reduce the number of polygons!rwa_adm_simple <- gSimplify(rwa_adm, tol=0.003) # simplifies the polygon # Resulting layer from gSimplify has no variables, so need to add a variable back in rwa_adm_simple$OBJECTID <- rwa_adm$OBJECTID
# 1. Merge data together -------------------------------------------------------
rwa_adm_simple <- merge(rwa_adm_simple, rwa_veg, by="OBJECTID")
# 2. Make dataframe interpretable by ggplot ------------------------------------
rwa_adm_df <- tidy(rwa_adm_simple, region="OBJECTID")
rwa_adm_df <- merge(rwa_adm_df, rwa_adm_simple, by.x="id", by.y="OBJECTID")
# 3. Make Plot -----------------------------------------------------------------
veg_map <- ggplot() +
geom_polygon(data=rwa_adm_df, aes(x=long, y=lat, group=group),
fill="white", color="black",size=.1) +
geom_polygon(data=rwa_adm_df[!is.na(rwa_adm_df$vegetation),],
aes(x=long, y=lat, group=group, fill=vegetation),
color="black",size=.1) +
theme_void() +
scale_fill_distiller(palette="RdYlGn", direction=1) +
labs(fill="Vegetation", title="Vegetation in Rwanda") +
theme(plot.title = element_text(hjust = 0.5, face="bold"))
veg_map
Satellite imagery often comes in multiple bands. Each band captures a different part of the electromagnetic spectrum. Different combinations of bands are used to create spectral indices which are used to highlight different land cover elements, such as vegetation, built-up area, or burnt areas.
\[NDVI = \frac{NIR-Red}{NIR+Red}\]
Landsat: Captures images across the earth every 16 days at a 30 meter resolution across multiple spectral bands. Landsat began collecting images in 1972 and continues to the present.
Sentinel: Captures images across the earth every 5 days (at the equator) at a 10 meter resolution across multiple sepctral bands. Began in 2014 and continues to the present.
MODIS: Captures daily to bi-weekly imagery ranging from a 250m to 500m resolution.
VIIRS: Monthly nighttime lights at a 750m resolution from April 2012 to present.
ndvi_2012_a <- raster(file.path(finalData, "Tiffs", "ndvi_2012_0203.tif")) plot(ndvi_2012_a) plot(ag_fields,add=T)
ndvi_2012_a <- crop(ndvi_2012_a, ag_fields) plot(ndvi_2012_a) plot(ag_fields, add=T)
ndvi_2012_a <- mask(ndvi_2012_a, ag_fields) plot(ndvi_2012_a)
You can also plot rasters using ggplot. Here, we need to convert the raster into a dataframe with latitude, longitude and the values. Then, we use geom_tile.
ndvi_2012_a_df <- as(ndvi_2012_a, "SpatialPixelsDataFrame") ndvi_2012_a_df <- as.data.frame(ndvi_2012_a_df) ggplot() + geom_tile(data=ndvi_2012_a_df, aes(x=x, y=y, fill=ndvi_2012_0203),alpha=1) + scale_fill_gradient(low="red",high="green") + coord_quickmap() + theme_void() + labs(fill="NDVI", title="NDVI")
We can also plot the raster using other packages that don't require first converting the object to a dataframe. Here, we use rastervis.
library(rasterVis)
cols <- colorRampPalette(brewer.pal(9,"RdYlGn"))
levelplot(ndvi_2012_a,
col.regions = cols)
levelplot(ndvi_2012_a,
col.regions = cols,
margin=F, # No histograms on side
par.settings=
list(axis.line=list(col='transparent')), # No box around image
scales=list(draw=FALSE), # No latitude and longitude markers
xlab="",
ylab="",
main="NDVI 2012 (Season A)",
colorkey=list(space="right"), # legend location. can change to bottom,
# left and top. If you don't want a
# legend, use: colorkey = F
maxpixels=1e6 # sometimes blends pixels together
) +
layer(sp.polygons(ag_fields,
col="black", # border color
lwd=1, # line width
alpha=1))
Note that NDVI goes from -1 to 1; MODIS scales this value by 10000 so we'll scale this back to -1 to 1 by multiplying values by 0.0001. There are two ways to do arithmetic on rasters:
raster_new <- raster * 0.0001
OR
raster_new <- calc(raster, fun = function(x) x * 0.0001)
The second way looks ugly, but it is more efficient.
library(rasterVis)
library(RColorBrewer)
ndvi_2012_a <- raster(file.path(finalData, "Tiffs", "ndvi_2012_0203.tif")) %>%
crop(ag_fields) %>%
mask(ag_fields) %>%
calc(fun=function(x) x*0.0001)
ndvi_2012_b <- raster(file.path(finalData, "Tiffs", "ndvi_2012_0910.tif")) %>%
crop(ag_fields) %>%
mask(ag_fields) %>%
calc(fun=function(x) x*0.0001)
ndvi_2016_a <- raster(file.path(finalData, "Tiffs", "ndvi_2016_0203.tif")) %>%
crop(ag_fields) %>%
mask(ag_fields) %>%
calc(fun=function(x) x*0.0001)
ndvi_2016_b <- raster(file.path(finalData, "Tiffs", "ndvi_2016_0910.tif")) %>%
crop(ag_fields) %>%
mask(ag_fields) %>%
calc(fun=function(x) x*0.0001)
ndvi_stack <- stack(ndvi_2012_a, ndvi_2012_b,
ndvi_2016_a, ndvi_2016_b)
cols <- colorRampPalette(brewer.pal(9,"RdYlGn"))
levelplot(ndvi_stack,
main="NDVI Across Seasons",
col.regions=cols,
layout=c(2, 2))
pal <- colorNumeric(c("red3", "yellow", "palegreen4"), values(ndvi_stack),
na.color = "transparent")
NDVI_interactive <- leaflet() %>%
addProviderTiles("OpenStreetMap") %>%
addRasterImage(ndvi_2012_a, colors = pal, opacity = 1, group="2012 A") %>%
addRasterImage(ndvi_2012_b, colors = pal, opacity = 1, group="2012 B") %>%
addRasterImage(ndvi_2016_a, colors = pal, opacity = 1, group="2016 A") %>%
addRasterImage(ndvi_2016_b, colors = pal, opacity = 1, group="2016 B") %>%
addLegend(pal = pal, values = values(ndvi_stack),
title = "NDVI") %>%
addLayersControl(
overlayGroups = c("2012 A","2012 B","2016 A","2016 B"),
options = layersControlOptions(collapsed = FALSE))
NDVI_interactive
Let's say we want to determine average NDVI values for each plot or the NDVI value under each household location. To do this, we use extract from the raster package.
# Average NDVI Value per Plot ag_fields$NDVI_2012_a <- extract(ndvi_2012_a, ag_fields, fun=mean) # NDVI Value per HH Location hh_data_sdf$NDVI_2012_a <- extract(ndvi_2012_a, hh_data_sdf)
Sometimes this process can be slow. To (dramatically) increase the speed of this task, use the velox package. This is built to work with spatial polygons (not points).
velox(raster) converts the raster object to a velox raster object. From the raster object you can use the extract method, where the inputs are the shapefile and summarize function.library(velox) ag_fields$ndvi_2012_a_velox <- velox(ndvi_2012_a)$extract(ag_fields, fun=mean)
subset(ag_fields@data, select=c(site, NDVI_2012_a, ndvi_2012_a_velox))
## site NDVI_2012_a ndvi_2012_a_velox ## 1 Kayonza15 0.4951754 0.4951754 ## 2 Kayonza4 0.5363880 0.5363880 ## 3 Rwamagana2 0.5958732 0.5958732 ## 4 Rwamagana33 0.5472275 0.5472275 ## 5 Rwamagana34 0.5667049 0.5667049 ## 6 Rwamagana35 0.5793644 0.5793644
Geographic Coordinate System
The most common geographic coordinate system is the World Geodetic System (WGS84).
"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
Projected Coordinate Systems
Below are projections for equal area and equal distance. Replace [LATITUDE VALUE] and [LONGITUDE VALUE] with the center of the area that you're working with.
Equal Area: Lambert Azimuthal Equal Area
"+proj=laea +lat_0=[LATITUDE VALUE] +lon_0=[LONGITUDE VALUE]"
Equal Distance: Azimuthal equidistant projection
"+proj=aeqd +lat_0=[LATITUDE VALUE] +lon_0=[LONGITUDE VALUE]"
A common task in GIS is to calculate distances. Here, let's calculate the distance of each household location to the nearest city (using a dataset of Rwanda's 3 largest cities). This same syntax can be used to calculate the shortest distance to polygons (e.g., a lake) or polylines (e.g., roads).
Import city-level dataset
rwa_cities <- read.csv(file.path(finalData, "rwa_cities.csv"))
coordinates(rwa_cities) <- ~lon+lat
crs(rwa_cities) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
Re-project to equal distant projection
# Here, I used Kigali as the center eq_dist_projection <- "+proj=aeqd +lat_0=-1.943889 +lon_0=30.059444" rwa_cities_ed <- spTransform(rwa_cities, CRS(eq_dist_projection)) hh_data_sdf_ed <- spTransform(hh_data_sdf, CRS(eq_dist_projection))
Creates distance matrix
gDistance(rwa_cities_ed, hh_data_sdf_ed, byid=T) %>% head
## 1 2 3 ## 1 39419.69 101785.29 74845.93 ## 2 41335.03 100056.47 75709.40 ## 3 41433.12 101555.78 76291.83 ## 4 40642.18 101794.90 75775.02 ## 5 40443.28 103553.65 76166.05 ## 6 39694.79 98492.55 73943.48
Calculate shortest distance to any city.
gDistance_cities_i <- function(i) gDistance(hh_data_sdf_ed[i,], rwa_cities_ed, byid=F) hh_data_sdf_ed$dist_city <- lapply(1:nrow(hh_data_sdf_ed), gDistance_cities_i) %>% unlist
There are also ways to calculate distances by using a geographic coordinate system (ie, keeping the world a sphere) and taking into account the curvator of the earth. Haversine (good for short distances – think the distance of a short haul flight) and vincenty (slower but better for longer distances) are common formulas to calculate these distances.
The below code uses the haversine method to calculate the distance between all the household locations and the first Rwandan city. The output is in meters.
library(geosphere) distHaversine(hh_data_sdf, rwa_cities[1,])
We can also calculate the distance between a point and a line using the dist2Line function from the geosphere package. Below shows an example between an arbitrary point and line shapefile using the haversine formula.
dist2Line(points, line, distfun=distHaversine)